#include "Polynomial.h"

#include <cmath>
#include <climits>
#include <cfloat>

#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define ACCURACY 1e-12

using namespace XYF;

// 显式构造函数
Polynomial::Polynomial(int o) : m_coefficients(std::vector<Real>(o, 0))
{
}

// 初始化列表构造
Polynomial::Polynomial(std::initializer_list<Real> l)
{
    for (auto it = l.begin(); it != l.end();)
        m_coefficients.push_back(*it++);
}

// 使用向量构造
Polynomial::Polynomial(const std::vector<Real> &v) : m_coefficients(v)
{
}

// 拷贝构造
Polynomial::Polynomial(const Polynomial &poly) : m_coefficients(poly.m_coefficients)
{
}

// 多项式赋值
void Polynomial::operator=(const Polynomial &poly)
{
    m_coefficients = poly.m_coefficients;
}

// 加减法
Polynomial Polynomial::operator+(const Polynomial &poly) const
{
    int index = 0;
    int O1 = order();
    int O2 = poly.order();

    // 取最大和最小的阶数
    int Maxsize = MAX(O1, O2);
    int Minsize = MIN(O1, O2);

    // 新的多项式阶数是最大的那个
    Polynomial res(Maxsize);

    // 先计算两者都有的次数的项
    while (index < Minsize)
    {
        res.m_coefficients[index] = m_coefficients[index] + poly.m_coefficients[index];
        index++;
    }

    // 添加补充的次数的项
    while (index < O1)
    {
        res.m_coefficients[index] = m_coefficients[index];
        index++;
    }
    while (index < O2)
    {
        res.m_coefficients[index] = res.m_coefficients[index] + poly.m_coefficients[index];
        index++;
    }

    return res;
}
Polynomial Polynomial::operator-(const Polynomial &poly) const
{
    int index = 0;
    int O1 = order();
    int O2 = poly.order();

    // 取最大和最小的阶数
    int Maxsize = MAX(O1, O2);
    int Minsize = MIN(O1, O2);

    // 新的多项式阶数是最大的那个
    Polynomial res(Maxsize);

    // 先计算两者都有的次数的项
    while (index < Minsize)
    {
        res.m_coefficients[index] = m_coefficients[index] - poly.m_coefficients[index];
        index++;
    }

    // 添加补充的次数的项
    while (index < O1)
    {
        res.m_coefficients[index] = m_coefficients[index];
        index++;
    }
    while (index < O2)
    {
        res.m_coefficients[index] = res.m_coefficients[index] - poly.m_coefficients[index];
        index++;
    }

    return res;
}

// 乘常数
Polynomial Polynomial::operator*(Real k) const
{
    int O = order();
    Polynomial res(O);
    for (int i = 0; i < O; i++)
        res.m_coefficients[i] = m_coefficients[i] * k;

    return res;
}
Polynomial Polynomial::operator/(Real k) const
{
    assert(fabs(k) > ACCURACY);

    int O = order();
    Polynomial res(O);
    for (int i = 0; i < O; i++)
        res.m_coefficients[i] = m_coefficients[i] / k;

    return res;
}

// 多项式相乘
Polynomial Polynomial::operator*(const Polynomial &poly) const
{
    const int O1 = order();
    const int O2 = poly.order();
    Polynomial res(O1 + O2 - 1);
    for (int i = 0; i < O1; i++)
        for (int j = 0; j < O2; j++)
            res.m_coefficients[i + j] += m_coefficients[i] * poly.m_coefficients[j];

    return res;
}

// 带余除法，获得商多项式
Polynomial Polynomial::operator/(const Polynomial &poly) const
{
    int O1 = order();
    int O2 = poly.order();

    // 阶数不够，返回 0 多项式
    if (O1 < O2)
        return Polynomial(0);

    Polynomial res(O1 - O2 + 1);

    // 最高阶系数
    res.m_coefficients[O1 - O2] = m_coefficients[O1 - 1] / poly.m_coefficients[O2 - 1];

    // 获得余式
    Polynomial tmp = *this - res * poly;

    // 减去最高项后需要降阶，然后递归计算商
    tmp.reduce();
    tmp = tmp / poly;

    // 商式赋值
    for (int i = 0; i < tmp.order(); ++i)
        res.m_coefficients[i] = tmp.m_coefficients[i];

    return res;
}

// 多项式求值
Real Polynomial::operator()(Real x) const
{
    // 先排除零多项式
    int O = order();
    if (O == 0)
        return 0;

    // 获得最高次系数，使用累乘方法提高效率
    Real y = m_coefficients.back();
    for (int i = O - 2; i >= 0; i--)
        y = y * x + m_coefficients[i];

    return y;
}

// 推入高阶项
void Polynomial::push(Real k)
{
    m_coefficients.push_back(k);
}

// 多项式降阶，将系数为 0 的高阶项去除
void Polynomial::reduce()
{
    // 只要最后一项非常小，就弹出最后一个元素，除非没有元素
    while (order() > 0 && fabs(m_coefficients.back()) < ACCURACY)
        m_coefficients.pop_back();
}

// 定积分
Real Polynomial::intergral(Real a, Real b) const
{
    int O = order();
    Real I = 0;
    for (int i = 0; i < O; i++)
        I += m_coefficients[i] * (pow(b, i + 1) - pow(a, i + 1)) / (i + 1);

    return I;
}

// 不定积分
Polynomial Polynomial::intergral() const
{
    Polynomial poly;
    poly.push(0);
    for (int i = 0; i < order(); i++)
        poly.push(coeff(i) / (i + 1));

    return poly;
}

// 导函数
Polynomial Polynomial::derivative() const
{
    Polynomial poly;
    for (int i = 1; i < order(); i++)
        poly.push(coeff(i) * i);

    return poly;
}

// 获取多项式的阶
int Polynomial::order() const
{
    return m_coefficients.size();
}

// 获得对应次数的系数
Real Polynomial::coeff(int n) const
{
    assert(n < m_coefficients.size());
    return m_coefficients[n];
}

// 多项式的拉普拉斯阵
Matrix Polynomial::Laplace() const
{
    int n = order() - 1;
    Matrix A(n, n);
    for (int i = 0; i < n; i++)
    {
        A(0, n - i - 1) = -m_coefficients[i];
        if (i < n - 1)
            A(i + 1, i) = 1;
    }
    return A;
}

// 没有区间端点，只有一个多项式
PiecewisePolynomial::PiecewisePolynomial(Polynomial poly) : m_poly({poly})
{
}

// 拷贝构造
PiecewisePolynomial::PiecewisePolynomial(const PiecewisePolynomial &ppoly) : m_interval(ppoly.m_interval), m_poly(ppoly.m_poly)
{
}

// 只设置区间，每段区间上是零
PiecewisePolynomial::PiecewisePolynomial(const std::deque<Real> &interval) : m_interval(interval), m_poly(std::deque<Polynomial>(interval.size() + 1))
{
    for (int i = 0; i < m_poly.size(); i++)
        m_poly[i] = {0};
}

// 设置区间和多项式
PiecewisePolynomial::PiecewisePolynomial(const std::deque<Real> &interval, const std::deque<Polynomial> &poly) : m_interval(interval), m_poly(poly)
{
    // 检查赋值是否有效
    assert(interval.size() + 1 == poly.size());
}

// 多项式赋值
void PiecewisePolynomial::operator=(const PiecewisePolynomial &poly)
{
    m_interval = poly.m_interval;
    m_poly = poly.m_poly;
}

// 多项式加减法
PiecewisePolynomial PiecewisePolynomial::operator+(const PiecewisePolynomial &poly) const
{
    // 必须具有相同的区间端点
    assert(m_interval.size() == poly.m_interval.size());
    for (int i = 0; i < m_interval.size(); ++i)
        assert(fabs(m_interval[i] - poly.m_interval[i]) < ACCURACY);

    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] + poly.m_poly[i];

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator-(const PiecewisePolynomial &poly) const
{
    // 必须具有相同的区间端点
    assert(m_interval.size() == poly.m_interval.size());
    for (int i = 0; i < m_interval.size(); ++i)
        assert(fabs(m_interval[i] - poly.m_interval[i]) < ACCURACY);

    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] - poly.m_poly[i];

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator+(const Polynomial &poly) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] + poly;

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator-(const Polynomial &poly) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] - poly;

    return res;
}

// 多项式乘常数
PiecewisePolynomial PiecewisePolynomial::operator*(Real k) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] * k;

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator/(Real k) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] / k;

    return res;
}

// 多项式相乘，分段相乘
PiecewisePolynomial PiecewisePolynomial::operator*(const PiecewisePolynomial &poly) const
{
    // 否则必须具有相同的区间端点
    assert(m_interval.size() == poly.m_interval.size());
    for (int i = 0; i < m_interval.size(); ++i)
        assert(fabs(m_interval[i] - poly.m_interval[i]) < ACCURACY);

    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] * poly.m_poly[i];

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator*(const Polynomial &poly) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] * poly;

    return res;
}

// 多项式除法
PiecewisePolynomial PiecewisePolynomial::operator/(const Polynomial &poly) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] / poly;

    return res;
}

// 多项式求值
Real PiecewisePolynomial::operator()(Real x) const
{
    return poly(x)(x);
}

// 在右端添加端点，增加一个零多项式
void PiecewisePolynomial::push_back(Real t)
{
    // 如果区间端点不是空，就检查是否大于最大的端点，如果重合或小于，就什么也不做
    if (!m_interval.empty() && t < m_interval.back() + ACCURACY)
        return;

    m_interval.push_back(t);
    m_poly.push_back(Polynomial());
}

// 在左端添加端点，增加一个零多项式
void PiecewisePolynomial::push_front(Real t)
{
    if (!m_interval.empty() && t > m_interval.front() - ACCURACY)
        return;

    m_interval.push_front(t);
    m_poly.push_front(Polynomial());
}

// 获得指定索引的多项式
Polynomial PiecewisePolynomial::operator[](int i) const
{
    assert(i < m_poly.size());
    return m_poly[i];
}

// 获得指定索引的多项式引用
Polynomial &PiecewisePolynomial::operator[](int i)
{
    assert(i < m_poly.size());
    return m_poly[i];
}

// 获得传入点对应的多项式
Polynomial PiecewisePolynomial::poly(Real t) const
{
    // 先排除超出区间范围的部分
    if (t < m_interval.front())
        return m_poly.front();
    if (t >= m_interval.back())
        return m_poly.back();

    // 找到 x 所属的区间范围的多项式，使用二分查找
    int a = 0, b = m_interval.size();
    int c = (a + b) / 2;
    // 注意左闭右开
    while (a + 1 < b)
    {
        if (t >= m_interval[c])
            a = c;
        else
            b = c;
        c = (a + b) / 2;
    }

    // 注意要排除最左边的那个多项式，所以加 1
    return m_poly[c + 1];
}

// 获得端点数
int PiecewisePolynomial::size() const
{
    return m_interval.size();
}

// 获得指定索引的阶数
int PiecewisePolynomial::order(int i) const
{
    assert(i < m_poly.size());
    return m_poly[i].order();
}

// 不定积分
PiecewisePolynomial PiecewisePolynomial::intergral() const
{
    PiecewisePolynomial ppoly;
    ppoly.m_interval = m_interval;

    // 注意默认构造的分段多项式有一个零多项式，需要清除
    ppoly.m_poly.clear();
    for (int i = 0; i < m_poly.size(); i++)
        ppoly.m_poly.push_back(m_poly[i].intergral());

    return ppoly;
}

// 导函数
PiecewisePolynomial PiecewisePolynomial::derivative() const
{
    PiecewisePolynomial ppoly;
    ppoly.m_interval = m_interval;

    // 注意默认构造的分段多项式有一个零多项式，需要清除
    ppoly.m_poly.clear();
    for (int i = 0; i < m_poly.size(); i++)
        ppoly.m_poly.push_back(m_poly[i].derivative());

    return ppoly;
}

// 检查端点处的连续性
int PiecewisePolynomial::continuity(int i) const
{
    assert(i < m_interval.size());

    // 左右两个多项式
    Polynomial left = m_poly[i];
    Polynomial right = m_poly[i + 1];

    int min = MIN(left.order(), right.order());
    Real t = m_interval[i];

    // 标记连续性，-1 表示不连续
    int c = -1;
    for (int j = 0; j < min; j++)
    {
        if (fabs(left(t) - right(t)) < ACCURACY)
            c++;
        left = left.derivative();
        right = right.derivative();
    }

    return c;
}

// 整体连续性
int PiecewisePolynomial::continuity() const
{
    // 不允许检查纯多项式的连续性，其连续性直接是无穷大
    assert(m_interval.size() > 0);

    int c = continuity(0);
    for (int i = 1; i < m_interval.size(); i++)
        c = MIN(c, continuity(i));
    return c;
}

// 合并两个分段多项式
PiecewisePolynomial XYF::merge(const PiecewisePolynomial &p1, const PiecewisePolynomial &p2)
{
    // 两个都是纯多项式，没有端点，直接相加
    if (p1.size() == 0 && p2.size() == 0)
        return p1 + p2;

    // 首先合并区间端点
    PiecewisePolynomial res;
    int n1 = p1.size();
    int n2 = p2.size();

    // 储存所有端点，然后排序
    Vector x(n1 + n2);
    for (int i = 0; i < n1; ++i)
        x[i] = p1.m_interval[i];
    for (int i = 0; i < n2; ++i)
        x[i + n1] = p2.m_interval[i];
    x.sort();

    // 推入这些端点
    for (int i = 0; i < n1 + n2; ++i)
        res.push_back(x[i]);

    // 开始计算加法，我们保存每个区间的中点，左右添加两个点
    std::deque<Real> mid;
    for (int i = 0; i < res.size() - 1; i++)
        mid.push_back((res.m_interval[i + 1] + res.m_interval[i]) / 2);
    mid.push_front(res.m_interval.front() - 1);
    mid.push_back(res.m_interval.back() + 1);

    // 获取对应的多项式相加
    for (int i = 0; i < mid.size(); i++)
        res.m_poly[i] = p1.poly(mid[i]) + p2.poly(mid[i]);

    return res;
}

// 阶乘
Real XYF::factorial(int n)
{
    Real i = 1;
    for (int k = 0; k < n; k++)
        i *= k + 1;

    return i;
}

// 组合数
Real XYF::combine(int n, int i)
{
    return factorial(n) / (factorial(i) * factorial(n - i));
}

// 多项式整数幂
Polynomial XYF::power(const Polynomial &poly, int n)
{
    Polynomial res({1});
    for (int i = 0; i < n; i++)
        res = res * poly;
    return res;
}

// 产生节点多项式
Polynomial XYF::nodes(const Vector &c)
{
    Polynomial poly({1});
    for (int i = 0; i < c.size(); i++)
        poly = poly * Polynomial({-c[i], 1});

    return poly;
}

// 已知节点 c0 c1 ... cn ，获得 cj 对应的初等多项式
Polynomial XYF::Lagrange(const Vector &c, int j)
{
    Polynomial poly({1});
    for (int k = 0; k < c.size(); k++)
        if (k != j)
            poly = poly * Polynomial({-c[k], 1}) / (c[j] - c[k]);

    return poly;
}

// Gauss-Legendre 多项式
Polynomial XYF::Legendre(int s)
{
    Polynomial F0({1});
    if (s > 0)
    {
        Polynomial x({0, 1});
        Polynomial F1 = x;
        Polynomial F2 = x;
        for (int i = 1; i < s; i++)
        {
            F2 = x * F1 * (2 * i + 1) / (i + 1) - F0 * i / (i + 1);
            F0 = F1;
            F1 = F2;
        }
        F0 = F2;
    }
    return F0;
}

// (0,1) 上的首一平移 Gauss-Legendre 多项式
Polynomial XYF::MonicLegendre(int s)
{
    Polynomial poly;
    for (int j = 0; j <= s; j++)
    {
        std::vector<Real> v(j + 1);
        v[j] = ((s - j) % 2 == 0 ? 1 : -1) * combine(s, j) * combine(s + j, j);
        poly = poly + Polynomial(v);
    }
    poly = poly * (pow(factorial(s), 2) / factorial(2 * s));
    return poly;
}

// 获得 Bernstein 基函数 Cni * t^i * (1-t)^(n-i)
Polynomial XYF::Bernstein(int n, int i)
{
    assert(n >= i);
    return power(Polynomial({0, 1}), i) * power(Polynomial({1, -1}), n - i) * combine(n, i);
}

// 获得 B 样条基函数
PiecewisePolynomial XYF::BasicBSpline(const std::deque<Real> &knot)
{
    // 至少要两个节点才能计算
    assert(knot.size() > 1);

    // 递归到最底层，恰好有两个节点
    if (knot.size() == 2)
    {
        PiecewisePolynomial Bik(knot);
        // 如果不重合，设置节点之间取常数 1
        if (knot[0] != knot[1])
            Bik[1] = Polynomial({1});
        return Bik;
    }

    // 减少一个节点得到第一阶的两个样条
    std::deque<Real> knot1 = knot;
    std::deque<Real> knot2 = knot;
    knot1.pop_back();
    knot2.pop_front();

    // 计算低阶样条
    PiecewisePolynomial Bik_1 = BasicBSpline(knot1);
    PiecewisePolynomial Bi_1k_1 = BasicBSpline(knot2);

    // 递推公式 B_i,k = (t-t_i)/(t_i+k-1-t_i)B_i,k-1 + (t_i+k - t)/(t_i+k - t_i+1)B_i+1,k-1
    Real t_i = knot[0], t_i1 = knot[1];
    Real t_ik_1 = knot[knot.size() - 2], t_ik = knot[knot.size() - 1];

    // 如果分母为零，则规定系数为零
    Real coeff1 = (fabs(t_ik_1 - t_i) < ACCURACY) ? 0 : (1 / (t_ik_1 - t_i));
    Real coeff2 = (fabs(t_ik - t_i1) < ACCURACY) ? 0 : (1 / (t_ik - t_i1));

    // 最后合并两个多项式
    return merge(Bik_1 * Polynomial({-t_i, 1}) * coeff1, Bi_1k_1 * Polynomial({t_ik, -1}) * coeff2);
}

// 计算部分差分表，宽度为 width
Table XYF::calculateTable(const std::vector<Real> &x, const std::vector<Vector> &f, int width)
{
    // 存放差分表
    Table table;

    // 计算差分表
    int index = 0;
    for (int i = 0; i < x.size();)
    {
        // 获得当前点的导数规模，也就是节点的重数
        int n = f[index].size();

        // 计算 n 重节点
        for (int j = 0; j < n; j++)
        {
            // 推入点值，然后递推计算
            table.push_back({f[index][0]});

            // 计算 i 行差分，不超过给定宽度
            for (int k = 1; k < MIN(width, i + 1); k++)
            {
                // 相邻的节点
                Real left = x[i - k];
                Real right = x[i];

                // 非重节点，直接差分计算
                if (left != right)
                    table[i].push_back((table[i][k - 1] - table[i - 1][k - 1]) / (right - left));
                else
                    // 如果遇到重节点，就推入高阶导数
                    table[i].push_back(1.0 / factorial(k) * f[index][k]);
            }

            // 换行
            i++;
        }

        // 记录不重复节点的索引
        index++;
    }

    return table;
}

// 多项式插值
Polynomial XYF::interpolate(const Vector &x, const Vector &f)
{
    int n = f.size();
    Table table; // 差分表

    // 计算差分表
    for (int i = 0; i < n; i++)
    {
        table.push_back({f[i]});
        for (int j = 1; j <= i; j++)
            table[i].push_back((table[i][j - 1] - table[i - 1][j - 1]) / (x[i] - x[i - j]));
    }

    // 由差分表求取插值多项式
    Polynomial poly(n);
    for (int i = 0; i < table.size(); i++)
    {
        Polynomial tmp({table[i][i]});

        // 多项式连乘
        for (int j = i - 1; j >= 0; j--)
            tmp = tmp * Polynomial({-x[j], 1});

        poly = poly + tmp;
    }

    return poly;
}

// NevilleAitken 直接计算插值点的值
Real XYF::NevilleAitken(const Vector &x, const Vector &f, Real t)
{
    // 存放差分表
    int N = x.size();
    Table table(N, std::vector<Real>(N));

    // 从下向上计算
    for (int i = N - 1; i >= 0; i--)
    {
        // 第一列直接赋值
        table[i][0] = f[i];

        // 从左向右计算
        for (int j = 1; i + j < N; j++)
            table[i][j] = ((t - x[i]) * table[i + 1][j - 1] -
                           (t - x[i + j]) * table[i][j - 1]) /
                          (x[i + j] - x[i]);
    }

    return table[0][N - 1];
}

// Hermite 牛顿插值
Polynomial XYF::NewtonInterp(const Vector &x0, const std::vector<Vector> &f)
{
    // 首先要把 x 进行扩充，添加重数
    std::vector<Real> x;
    for (int i = 0; i < f.size(); i++)
    {
        // 获得当前点的导数规模，也就是节点的重数
        int n = f[i].size();

        for (int j = 0; j < n; j++)
            x.push_back(x0[i]);
    }

    // 计算差分表
    Table table = calculateTable(x, f, x.size());

    // 由差分表求取插值多项式
    Polynomial poly;
    for (int i = 0; i < table.size(); i++)
    {
        Polynomial tmp({table[i][i]});

        // 多项式连乘
        for (int j = i - 1; j >= 0; j--)
            tmp = tmp * Polynomial({-x[j], 1});

        poly = poly + tmp;
    }

    return poly;
}

// 普通 3 次样条插值
PiecewisePolynomial XYF::ppFormSpline(const Vector &x, const std::vector<Vector> &f, BCType type)
{
    // 差分表
    int N = x.size();

    // 将插值条件作为区间初始化
    PiecewisePolynomial ppoly;
    for (int i = 0; i < N; ++i)
        ppoly.push_back(x[i]);

    // 首先要把 x 进行扩充，添加重数
    std::vector<Real> y;
    for (int i = 0; i < f.size(); i++)
    {
        // 获得当前点的导数规模，也就是节点的重数
        int n = f[i].size();

        for (int j = 0; j < n; j++)
            y.push_back(x[i]);
    }

    // 计算差分表，宽度不超过 3，因为是 3 次样条
    Table table = calculateTable(y, f, 3);

    // 初始化线性系统
    Matrix A(N, N);
    Vector B(N), lambda(N - 2), mu(N - 2);

    // 计算系数
    for (int i = 0; i < N - 2; i++)
    {
        mu[i] = (x[i + 1] - x[i]) / (x[i + 2] - x[i]);
        lambda[i] = (x[i + 2] - x[i + 1]) / (x[i + 2] - x[i]);
    }
    // 系数矩阵
    for (int i = 1; i < N - 1; i++)
    {
        A(i, i - 1) = mu[i - 1];
        A(i, i + 1) = lambda[i - 1];
        A(i, i) = 2;
    }

    // 根据样条类型修改系统
    switch (type)
    {
    // 完全样条
    case BCType::complete_spline:
    {
        A(0, 0) = 2, A(0, 1) = 1;
        A(N - 1, N - 2) = 1, A(N - 1, N - 1) = 2;

        for (int i = 0; i < N; i++)
            B[i] = 6 * table[i + 2][2];

        break;
    }
    // 无结样条
    case BCType::notaknot_spline:
    {
        // 利用三阶导相等
        A(0, 0) = x[2] - x[1];
        A(0, 1) = x[0] - x[2];
        A(0, 2) = x[1] - x[0];
        A(N - 1, N - 3) = x[N - 1] - x[N - 2];
        A(N - 1, N - 2) = x[N - 3] - x[N - 1];
        A(N - 1, N - 1) = x[N - 2] - x[N - 3];

        for (int i = 1; i < N - 1; i++)
            B[i] = 6 * table[i + 1][2];

        B[0] = B[N - 1] = 0;

        break;
    }
    // 周期样条
    case BCType::periodic_spline:
    {
        A(0, 0) = 1;
        A(0, N - 1) = -1;
        A(N - 1, 0) = 2 * (x[1] - x[0]);
        A(N - 1, 1) = x[1] - x[0];
        A(N - 1, N - 2) = x[N - 1] - x[N - 2];
        A(N - 1, N - 1) = 2 * (x[N - 1] - x[N - 2]);

        for (int i = 1; i < N - 1; i++)
            B[i] = 6 * table[i + 1][2];

        // 6(f[x_1, x_2] - f[x_N - x_{N-1}])
        B[0] = 0;
        B[N - 1] = 6 * (f[1][0] - f[0][0]) / (x[1] - x[0]) - (f[N - 1][0] - f[N - 2][0]) / (x[N - 1] - x[N - 2]);
        // 由于周期样条插值不需要重节点，因此节点和值一一对应

        break;
    }
    }

    // 求解二阶导的向量
    LinearSolver solver;
    Vector M = solver.solve(A, B, true);

    Real fPrime = 0;
    // 完全样条有一阶导信息
    if (type == BCType::complete_spline)
        fPrime = f[0][1];
    else
        // 根据二阶导向量补充一阶导的信息
        // f[x_1, x_1] = f[x_1, x_2] - (2M_1 + M_2) / 6 * (x_2 - x_1)
        fPrime = table[1][1] - (2 * M[0] + M[1]) / 6 * (x[1] - x[0]);

    // 利用导数条件推导多项式
    for (int i = 0; i < N - 1; i++)
    {
        // 泰勒展开的位置
        Real t = x[i];
        Polynomial c0({f[i][0]});
        Polynomial c1({-t, 1});

        // 利用泰勒展开式
        // s_i(x) = f_i + f[x_i, x_i](x - x_i) + M_i / 2 * (x - x_i)^2 + (M_{i+1} - M_i) / 6(x_{i+1} - x_i) * (x - x_i)^3
        Polynomial poly = c0 + c1 * fPrime + power(c1, 2) * M[i] / 2 + power(c1, 3) * (M[i + 1] - M[i]) / (6 * (x[i + 1] - x[i]));

        // 获得新的一阶导条件，同时确定此段区间上的多项式
        fPrime = poly.derivative()(x[i + 1]);
        ppoly[i + 1] = poly;
    }

    return ppoly;
}

// 2 次基数 B 样条插值
PiecewisePolynomial XYF::cardinalBSpline2(const Vector &x, const Vector &f)
{
    // 差分表
    int N = x.size();

    // 将插值条件作为区间初始化
    PiecewisePolynomial ppoly;
    for (int i = 0; i < N; ++i)
        ppoly.push_back(x[i]);

    // 注意一共有 N 个插值点，但两端都是整数点，需要排除
    // 插值点的个数是总数减一
    N = N - 1;
    Vector b(N - 1);
    Matrix M(N - 1, N - 1);

    // 系数多项式
    for (int i = 1; i < N - 2; i++)
    {
        M(i, i - 1) = M(i, i + 1) = 1;
        M(i, i) = 6;
    }
    M(0, 0) = M(N - 2, N - 2) = 5;
    M(0, 1) = M(N - 2, N - 3) = 1;

    // 目标向量
    for (int i = 1; i < N - 2; i++)
        b[i] = 8 * f[i + 1];

    b[0] = 8 * f[1] - 2 * f[0];
    b[N - 2] = 8 * f[N - 1] - 2 * f[N];

    // 获得样条系数
    LinearSolver solver;
    Vector a(N + 1);
    Vector tmp_a = solver.solve(M, b, true);
    for (int i = 0; i < N - 1; i++)
        a[i + 1] = tmp_a[i];

    a[0] = 2 * f[0] - a[1];
    a[N] = 2 * f[N] - a[N - 1];

    // 累计样条函数
    std::deque<Real> knot = {x[0] - 2, x[0] - 1, x[0], x[0] + 1};
    for (int i = 0; i < N + 1; i++)
    {
        PiecewisePolynomial BS = BasicBSpline(knot) * a[i];
        ppoly = merge(ppoly, BS);

        // 平移节点
        knot.pop_front();
        knot.push_back(knot.back() + 1);
    }

    return ppoly;
}

// 3 次基数 B 样条插值
PiecewisePolynomial XYF::cardinalBSpline3(const Vector &x, const std::vector<Vector> &f)
{
    // 差分表
    int N = x.size();

    // 将插值条件作为区间初始化
    PiecewisePolynomial ppoly;
    for (int i = 0; i < N; ++i)
        ppoly.push_back(x[i]);

    // 三次 B 样条，在整数点插值
    Vector B(N);
    Matrix M(N, N);

    // 系数多项式
    for (int i = 1; i < N - 1; i++)
    {
        M(i, i - 1) = M(i, i + 1) = 1;
        M(i, i) = 4;
    }
    M(0, 0) = M(N - 1, N - 1) = 4;
    M(0, 1) = M(N - 1, N - 2) = 2;

    // 目标向量
    for (int i = 1; i < N - 1; i++)
        B[i] = 6 * f[i][0];

    B[0] = 6 * f[0][0] + 2 * f[0][1];
    B[N - 1] = 6 * f[N - 1][0] - 2 * f[N - 1][1];

    // 求解获得样条系数
    LinearSolver solver;
    Vector a(N + 2);
    Vector tmp_a = solver.solve(M, B, true);

    for (int i = 0; i < N; i++)
        a[i + 1] = tmp_a[i];

    a[0] = a[2] - 2 * f[0][1];
    a[N + 1] = a[N - 1] + 2 * f[N - 1][1];

    // 累计样条函数
    std::deque<Real> knot = {x[0] - 3, x[0] - 2, x[0] - 1, x[0], x[0] + 1};
    for (int i = -1; i < N + 1; i++)
    {
        PiecewisePolynomial BS = BasicBSpline(knot) * a[i + 1];
        ppoly = merge(ppoly, BS);

        // 平移节点
        knot.pop_front();
        knot.push_back(knot.back() + 1);
    }

    return ppoly;
}
